1507.01946v2 [cond-mat.quant-gas] 15 Nov 2015 


e 
e 


arXiv 


Prethermal Floquet Steady States and Instabilities in the Periodically Driven, Weakly 
Interacting Bose-Hubbard Model 


Marin Bukov,! Sarang Gopalakrishnan,” Michael Knap,?’? and Eugene Demler? 


‘Department of Physics, Boston University, 590 Commonwealth Ave., Boston, MA 02215, USA 
? Department of Physics, Harvard University, 17 Oxford Street, Cambridge, MA 02138, USA 
*Department of Physics, Walter Schottky Institute, and Institute for Advanced Study, 
Technical University Munich, 85748 Garching, Germany 
(Dated: November 17, 2015) 


We explore prethermal Floquet steady states and instabilities of the weakly interacting two- 
dimensional Bose-Hubbard model subject to periodic driving. We develop a description of the 
nonequilibrium dynamics, at arbitrary drive strength and frequency, using a weak-coupling con- 
serving approximation. We establish the regimes in which conventional (zero-momentum) and un- 
conventional [(7,7)-momentum] condensates are stable on intermediate time scales. We find that 
condensate stability is enhanced by increasing the drive strength, because this decreases the band- 
width of quasiparticle excitations and thus impedes resonant absorption and heating. Our results 
are directly relevant to a number of current experiments with ultracold bosons. 


Periodically driven systems[I}{4] often exhibit ex- 
otic phenomena that are absent in their non-driven 
counterparts|[5}H7]. Classic examples include the Kapitza 
pendulum and the periodically kicked rotor. Recently, 
periodically modulating optical lattices has attracted in- 
terest as a way of controlling hopping processes[8}{14] in 
order to engineer gauge fields[15}{22|, topological band 
structures , and associated exotic states of matter. 
Such exotic states are known to exist in noninteracting 
systems and in certain mean-field models; the extent to 
which they survive in the presence of interactions is a 
central open question. It is believed, from the eigenstate 
thermalization hypothesis[30}32], that driven interacting 
systems will generically heat up to infinite temperature 
at sufficiently late times[33H40]. Nevertheless, in some 
parameter regimes these heating times will be paramet- 
rically slower than the system’s characteristic time scales. 
In that case, the system will rapidly approach a “prether- 
malized” Floquet steady state[39] [41}43], which governs 
the dynamics until the much later heating time scales. 

In the present work, we study these prethermal states 
in the weakly interacting, two-dimensional, periodically- 
driven Bose-Hubbard model (BHM). The regime we ex- 
plore is directly relevant to experiments{10} [IT] [14] [T7}+ 
[19] (27) [22] [26], in which weak interactions are present. 
We employ a self-consistent weak-coupling conserving ap- 
proximation (WCCA) which treats the coupled nonlinear 
dynamics of the condensate and the quasiparticle spec- 
trum while neglecting collisions between quasiparticles. 
This approximation is justified at weak coupling since 
nonlinearities are important at much shorter times than 
the collisional time scales. 

Within the WCCA, we find a phase diagram (Fig. |1) 
featuring at low drive frequency a regime in which the 
superfluid state is already unstable within Bogoliubov 
theory, owing to the resonant creation of quasiparticle 
pairs, and a regime (at high drive frequency) where the 
superfluid is stable. In the WCCA, there is a sharp 


phase transition between these; when effects beyond weak 
coupling are included, there will be a qualitative differ- 
ence in heating rates. Thus, in the “stable” regions of 
Fig. |1| the system initially reaches a prethermalized su- 
perfluid state—featuring a nonequilibrium quasiparticle 
distribution—and then eventually heats up. For strong 
driving, the prethermalized superfluid state is exotic, in- 
volving condensation at momentum ma = (7,7). The 
existence of this exotic phase in the high-frequency limit 
has previously been established {8} [9] [11]; we find that it 
persists for intermediate frequencies as well. 
Remarkably, we find that the stable phase is enhanced 
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FIG. 1: (Color online) Stability diagram of the driven BHM 
for U/Jo = 0.2. In the pink regions the condensate is unsta- 
ble as the drive parametrically excites pairs of quasiparticles. 
In contrast, in the blue regions the condensate is stable on 
intermediate time scales. In the grey shaded region around 
¢ & 2.405 the system is strongly correlated (see text). The 
symbols represent numerical WCCA results; the boundaries 
are given by the analytical expression Eq. (5). Points marked 
(a), (b), (c) correspond to the panels in Fig. 


for intermediate drive strengths, since the drive both cre- 
ates quasiparticle pairs when this is a resonant process, 
and decreases the effective hopping rate and thus the 
effective bandwidth of quasiparticle excitations. A key 
conclusion of our work is that, for weak interactions but 
general drive amplitude and frequency, the condensate 
becomes unstable when the drive frequency is parametri- 
cally resonant with the drive-renormalized time-averaged 
bandwidth. Therefore, parametric resonance occurs at 
lower frequencies when the drive strength is ramped up. 
Model.—We consider the Bose-Hubbard model on a 
square lattice in the presence of a circularly-polarised 
time-periodic force E(t) = A (cos M¢, sin Mt)” 


Map(t = Joy ibs ye [Sra 


The operator bi creates a boson on lattice site rj. The 
tunnelling and interaction strength are denoted by Jo and 
U, respectively. To achieve non-trivial dynamics in the 
high-frequency regime, we scale the driving amplitude 
linearly with the driving frequency A ~ 2.|6]; we define 
¢ = A/Q. We transform this Hamiltonian into a rotating 
frame (cf. supplementary material{75]), giving: 


zn . (2) 


Thus, in the rotating frame, the system experiences 
an effective time-dependent gauge potential A(t) = 
¢ (sin Nt, — cos Mt)” The time evolution of U(1)- 
invariant quantities (and thus the stability) remains the 
same in both frames[44]. 

Method. To study the driven system at arbitrary fre- 
quencies, we employ a self-consistent, weak-coupling con- 
serving approximation (WCCA). The WCCA involves 
deriving equations of motion from a two-particle irre- 
ducible effective action[45] within the nonequilibrium 
Schwinger-Keldysh formalism[46] [47], keeping only di- 
agrams to first order in U (see [75]). Unlike simple 
perturbation theory or Bogoliubov theory, the WCCA 
respects unitarity and conservation laws[48], and thus 
gives physically sensible results for all times; in partic- 
ular, it allows the exponential growth of unstable modes 
to be cut off by the resulting depletion of the conden- 
sate. While the WCCA is not guaranteed to yield a 
gapless excitation spectrum|48] [49], the low-frequency 
behavior of the spectrum is irrelevant for the phenom- 
ena discussed here. Our approach is equivalent to a fully 
self-consistent, time-dependent Hartree-Fock-Bogoliubov 
(HFB) approximation [49] {50}; our formulation, however, 
can more readily be extended to higher orders in U. 

The WCCA equations of motion[75] were solved nu- 
merically. For the results presented here, we prepared 
the system on a N, = 100 x 100 lattice in the ground 
state of Bogoliubov theory. We allow for a macroscopic 
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FIG. 2: (Color online). (a) Time evolution of the condensate 
fraction for 801 driving cycles, starting from a Bogoliubov ini- 
tial state localised at k = 0 for U/Jo = 0.2. (b) Decay rate to 
75% of the condensate curves for 2/Jo = 12 (boldface points 
in Fig. [1p. Error bars are set by the difference of the inverse 
times, determined by the first and last time the curve passes 
through 3/4 taking into account the oscillatory behaviour. 


population of the k = 7 mode to allow for a condensate 
at momentum 7. To study the nonequilibrium dynam- 
ics, we abruptly turn on the periodic drive and propagate 
the initial state for 801 driving cycles using Eqs. (15) and 
(16) off[75]. We checked that the results are insensitive 
to system size. 

Stability diagram.—The stability phase diagram is 
shown in Fig. Previous work has investigated 
the driven Bose-Hubbard model{5I}58) and related 
models[59}H67| using various approximation schemes; we 
go beyond these works by treating both the condensate 
and quasiparticle sectors, including the feedback between 
them. Thus, we are able to explore instabilities originat- 
ing in either sector on equal footing. 

We first discuss two analytically tractable limits, cor- 
responding to high-frequency driving (i.e., going along 
the x axis of Fig.|1) and to low-amplitude driving (i.e., 
going along the y axis). In the first case, the dynamics 
is approximately governed by an effective time-average 
Hamiltonian[5} [6]: 


aye = Java Od blb; “pS ay 5 (nj ~ 1). (3) 
J 


The periodic modulation renormalizes the hopping to 
Jave(¢) = JoJo(C), where Jo(¢) is the zeroth-order Bessel 
function of the first kind, which is a damped oscilla- 
tory function with the first zero at ¢ = 2.4, the sec- 
ond at ¢ = 5.5, etc. Thus, as ¢ is increased, the time- 
averaged hopping decreases, until the dispersion flattens 
at ¢ = 2.4. For ¢ > 2.4 the dispersion flips sign, and 
acquires a stable minimum at a = (7,7). Thus, in the 
high-frequency limit the condensate at 0 = (0,0) is stable 
when ¢ < 2.4, whereas the condensate at 7 is stable when 
2.4 S ¢ S 5.5. Moreover, for commensurate filling, the 
superfluid phase should transition into a Mott insulating 
state around ¢ = 2.4 determined by the phase boundary 
Jave(€)/U < 0.06. [68] [69] This transition regime, marked 
by the thin vertical strip in Fig. [1] is beyond the validity 
of the WCCA; our WCCA simulations in this regime give 
oscillatory behavior, see [75]. 

A second analytically tractable limit is that of weak 
driving, at arbitrary Q. The dominant effects can be 
inferred from linear stability analysis around the non- 
driven state. In terms of Bogoliubov quasiparticle op- 
erators yx, the system-drive coupling includes terms of 
the form etal at, involving the emission of pairs of 
quasiparticles from the condensate. The emission rate 
is related to the density of states of two-quasiparticle 
excitations at Q. Specifically, if the non-driven system 
has quasiparticle excitations at energy Ex, E_, such that 
Q = Ex, + E_,, absorption will occur and the system 
will be unstable. On the other hand, if Q > 2W, where 
W & 2zJp is the approximate bandwidth of Bogoliubov 
excitations, then absorption does not occur and the sys- 
tem is stable. 

Combining the insights from these two limits al- 
lows us to understand the entire stability phase dia- 
gram. The drive creates pairs of renormalized Bogoli- 
ubov quasiparticles, which have an effective bandwidth 
Wave © 2zJave(C). We define Wave = maxg[Fave(k)] — 
minx |Eave(k)] as the time-averaged Floquet-Bogoliubov 
bandwidth; in terms of this, the stability condition reads 


Q. > 2Wave & stable. (4) 


Equation (4) is consistent with our numerical re- 
sults (Fig. |1). This result is unexpected—since the 
time-averaged Hamiltonian is valid at infinite fre- 
quency whereas parametric resonance is a low-frequency 
phenomenon— but can be understood as follows. The 
hopping matrix element in the driven system can be ex- 
panded as J(t) ~ Jo >>, In(¢) exp(inQt). We absorb 
the time-independent n = 0 component in the unper- 
turbed Hamiltonian, and treat the n = 1 term, which 
oscillates at Q, perturbatively. The perturbation is small 
for U < Q, because the matrix element for creating 
two quasiparticles is proportional to both J(¢) [which 
need not be small] and U [which is assumed to be small]. 
We then use parametric instability analysis[75] with the 
renormalized dispersion, and conclude that an instability 


occurs when 2 = 2Waye. When Q/Jo >> 1, the critical 
driving frequency is given by 


2-(¢) = Aa/ 2Jeve(€)(Zave(C) + no). (5) 


Note that in the present case, resonant absorption occurs 
for drive strengths up to twice the single-particle band- 
width; by contrast, in noninteracting systems, no absorp- 
tion occurs for Q > Waye. The presence of absorption 
at frequencies exceeding the single-particle bandwidth is 
generic in interacting systems. 

Condensate evolution.— Figure [2 [panel (a)] shows the 
evolution of the condensate fraction in various regimes: 
in the parametrically unstable regime (solid blue line), 
the condensate slowly decays; in the stable regime 
(dashed red line), it saturates to a prethermalized value, 
which is generally lower than the Bogoliubov value (since 
|Jave(¢)| < |Jo|). The system enters a steady-state with 
constant in time evolution when measured stroboscopi- 
cally. When the initial condensate is at the band max- 
imum (dash-dotted black line), the condensate decays 
rapidly. Panel (b) shows the decay rate as a func- 
tion of drive amplitude in the parametrically unstable 
regime: note that the decay rate depends not only on 
drive strength ¢, but also on U and 2. Very close to 
the region ¢ ~ 2.405 (grey strip in Fig. [1), the WCCA 
gives strong oscillations of the particle density between 
the condensates at 0 and 7 (see [75]); however, as previ- 
ously noted, the WCCA is not reliable here. 

A natural further observable is the total energy of the 
system, which grows in the unstable phases and saturates 
in the stable phases (see [75}). 

(Quasi-)momentum distribution.—Fig. |3} plots snap- 
shots of the quasimomentum (i.e., lattice momentum) 
distribution; the time evolution of this quantity is shown 
in[75]. Specifically, the quantity plotted is nx, = (bt by) - 
nodk,0, i.e., the condensate peak is subtracted. The 
quasimomentum distribution can be directly accessed 
through band mapping followed by time-of-flight imag- 
ing. Moreover, as we are concerned with a single-band 
model, one can extract this distribution directly from 
time-of-flight imaging, by focusing on momenta within 
the first Brillouin zone. 

Figure |3] (a) shows the parametrically unstable case, 
where quasiparticles are strongly excited around the 
quasimomentum surface {k : Q = 2Faye(k)} matching 
the resonance condition. Within Bogoliubov theory, the 
(time-averaged) excitation intensity should be uniform 
along this surface. However, as the points along this 
surface are not symmetry-related, the nonlinearities in- 
cluded in the WCCA favor some points on the excitation 
surface, as seen in the intensity pattern in Fig. |3) (a). 

Figure [}| (b) shows the stable case. Here, by contrast 
with panel (a), the quasiparticle population remains low 
throughout the Brillouin zone. As expected from Bo- 
goliubov theory, bosonic modes satisfying Jaye(k) < U 
should have appreciable occupation in the steady state; 








-r/2 0 nl2 -n/2 0 m/2 T 
ka ky ky 


FIG. 3: (Color online). Snapshot of the momentum distribution nx = (b). bk) — nodx,o after 801 driving cycles starting from a 
Bogoliubov initial state localised at k = O for U/Jo = 0.2. Panel (a) is in the unstable regime where the condensate is depleted 
due to parametric resonance. The bosons are excited by the drive to the quasienergy surface 2 = 2Faye(k) (bright yellow-white 
circle around k = 7) where they occupy sharp peaks (white pixels). Panel (b) is in the regime where the condensate is stable on 
the pre-thermal time scales. In panel (c), the system is dynamically unstable due to the dispersion being inverted. The bright 
disc of excitations around k = 0 corresponds to dynamically unstable modes. The parameters are (a) 2./Jo = 10, ¢ = 0.8, (b) 


Q/Jo = 18, € = 2.2, and (c) Q/Jo = 20, € = 3.8. 


this region expands as the dispersion flattens. The intri- 
cate patterns in momentum space are due to the abrupt 
turn-on of the drive—which initializes the Floquet- 
Bogoliubov quasiparticle states out of equilibrium—and 
are absent when the drive is instead gradually ramped up. 
These patterns evolve nontrivially with time (see [75]). 
Finally, Fig. |3|(c) illustrates the case in which the ini- 
tial state is a condensate at k = 0, but the dispersion 
is inverted (¢ > 2.4) so that the only stable conden- 
sate is supported at k = zw. Thus the initial state is 
unstable regardless of 2. Let us consider the infinite- 
frequency limit; which amounts to a sudden quench of 
the single-particle dispersion. Computing the Bogoli- 


ubov spectrum around a condensate at k = O in an 
inverted dispersion, we find that modes with momenta 
near k = O acquire imaginary frequencies (and thus 


grow exponentially), whereas modes with large momenta 
are stable [81]. The unstable modes are determined by 
the condition Egye(k) + zJg < 2n9U, where €aye(k) is 
the single-particle Floquet dispersion (3). These modes 
are dynamically stabilized due to the nonlinear feedback 
of the self-consistent treatment{47]. Our numerical re- 
sults with the WCCA confirm this picture: the unstable 
modes at small quasimomenta acquire large populations, 
whereas the large-quasimomentum modes do not. This 
behavior is specific to the WCCA; in a real system it will 
correspond to intermediate-time dynamics t < Jo/U?. 
On longer times, collisions between quasiparticles should 
cause large occupation numbers across the Brillouin zone, 
see [75]. 

Discussion.—We briefly outline the validity of the 
WCCA in the three regimes of interest (for details 
see ). In the parametrically unstable regime, the sta- 
bility analysis suggests that unstable modes grow at the 
rate T ~ UngJoAi(C)/Wave, while the momentum arcs in 
Fig.|2|(a) decay at a Golden Rule rate ~ U?nonk/Wave- 


Hence, as long as Unx < JoZi(¢), the formation rate is 
greater than the decay rate and the WCCA is reliable. 
In the stable region, the condensate fraction no remains 
large, and the WCCA remains valid, until very late times, 
when resonant absorption involving m = 2/Waye quasi- 
particles becomes dominant. For large ©, this is a very 
high-order and therefore very slow process. Finally, in 
the dynamically unstable phase, the WCCA physics is 
valid up to times W,y./U? (the collisional time scale). 
Thus, at weak coupling, there is a parametrically large 
window between 1/\/Jave(C)U and Waye/U? where the 
WCCA description is correct. 


The main experimental prediction of this work—a 
parametric change in heating rates as a function of 
drive amplitude and frequency—can be measured in 
present-day experiments, which are naturally in the 
weak-coupling regime. For the experiment in Ref. [26] the 
parameters were chosen as U/Jp & 0.1, Q/Jo 20, and 
¢ & 0.6, which is within the regime we considered. For 
realistic experiments in optical lattices, the presence of 
higher bands can lead to instability even at high drive 
frequencies 2. In this case there are three regimes: (i) 
if 2 is less than twice the renormalized bandwidth W,,. 
of the lower band, the system is parametrically unsta- 
ble as discussed above; (ii) if Q is larger than 2Wave, 
smaller than the band gap to the upper band, and fur- 
thermore chosen such that any n-photon resonances to 
higher bands are suppressed, then the system is sta- 
ble within WCCA. (iii) if Q exceeds the band gap, the 
drive can mediate interband transitions, leading to in- 
stability again. For a square optical lattice with typical 
lattice potential Viatt = 10 Fy ecou, Evecoi! = h x 4 kHz, the 
bandwidth of the lowest band is Wo = 4Jp = h x 0.3 kHz 
[the time-averaged bandwidth W,,- is reduced by a fac- 
tor of Jo(¢)], and the gap to the second Bloch band is 
A = 4.57 Execoil) = h X 18.28 kHz. 


Although we focused on a square lattice, the argu- 
ments generalize to other lattices, such as the honey- 
comb lattice, in which topologically non-trivial states ex- 
ist. Note that topological gaps in mechanically shaken 
optical lattices scale as Q~'[23}25]. Hence, in order to 
engineer topological insulators with large gaps (and a 
large region of non-zero Berry curvature around them), 
it is desirable to go to lower frequencies. Our results im- 
pose a fundamental limit for weakly-interacting bosonic 
systems on how small the frequency can be, since for 
Q < 2Waye the system becomes unstable. More gener- 
ally, our results suggest that conserving approximations, 
whether controlled by weak coupling or some other pa- 
rameter as in large-N models|47 [71H74], are ways of ex- 
ploring dynamical phase transitions in models that are 
both interacting (unlike free-particle models) and _finite- 
dimensional (unlike the Kapitza pendulum). The critical 
properties of such transitions are a fruitful theme for fu- 
ture work. Although in practice such phase transitions 
will be smeared out by higher-order effects, the associated 
crossovers should still be experimentally observable. 
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This might seem counterintuitive, as the larger- 
momentum modes have “more negative” energies; note, 
however, that in the U —> 0 limit, all modes are stable 
as there are no decay processes. 





TRANSFORMATION TO THE ROTATING FRAME AND STABILITY ANALYSIS 


In this appendix, we begin by discussing the transformation of the driven BHM from the lab frame to the rotating 
frame, before doing a Bogoliubov stability analysis. The rotating frame is defined by the unitary transformation 


V(t) =exp | -i | | B(U}at! + A(0)| “xn, , {rot (t)) =V(t)|\dran(t)), bs 4 be AO, (6) 


where A(t) = ¢ (sin Qt, — cos at)”. This time-dependent change of basis is equivalent to the standard gauge trans- 


formation in electromagnetism E(t) ~ 0;A(t). 


Physically, the transformation trades the fast time-dependence of 


the quasimomentum for a time oscillating dispersion relation. The infinite-frequency limit is non-trivial when the 
amplitude of the gauge potential ¢ = A/Q remains finite. This rotation facilitates the analytic computation the 
time-averaged Floquet Hamiltonian. In fact, going to the rotating frame is equivalent to re-summing an infinite 
inverse-frequency subseries in the lab-frame]|6]. As a result, the effective hopping matrix element is a non-perturbative 
function of ¢. 


The rotating frame Hamiltonian reads 


H(t) 


(J) j 
= Ave ~ S- (Joe Ae) = Jess) bib, 
(13) 
where in the second equality we separated out the time average explicitly. The time-average Hamiltonian Haye is 
defined in Eq. (3) of the main text, and the effective hopping is Jaye = JoJo(¢). To avoid confusion in the notation, 
we note in passing that this effective Hamiltonian is not equivalent to the full Floquet Hamiltonian at any finite 
driving frequency, see eg. Refs. [5] [6]. 

Before we dive into the details of the parametric stability analysis for the driven Bose-Hubbard model, let us 
demonstrate how to derive the stability criterion with the help of the Rotating Wave Approximation (RWA) [a very 
similar method was used to study the parametric instability in periodically-driven Luttinger liquids [76]]. For this 
purpose, we choose the parametric oscillator with Hamiltonian H(t) = § (p? + w§a? + awé cos(Nt)x?). Writing the 
Hamiltonian using ladder operators x = 1/,/2wo(y! +7) and p = iy/wo/2(7! — 7), and dropping any (time-dependent) 
constants leads to 


H(t) = wo (1 + 5 0s or) aly + = cos Ot (yryt + hie.) . 


If we parametrise y(t) = u’(t)y(t = 0) — u’*(t)yT(t = 0), with u/(t = 0) = 1 and v(t = 0) = 0, we can write 


Heisenberg’s EOM as 
io u'\ (wo + $wo cos Nt Fw cos Qt ul 
dt \v J — —FuwocosNt — (wo + $wo cos Qt) v! 


= [woo? + W(t)] (ir) + Gancoser( ° - ee? (7) 


where the matrix W(t) = $wo cos(Qt)o* has zero time-average and o* is the Pauli matrix in Bogoliubov space. We 
now apply the transformation w(t) = e’?°'u’(t), o’(t) = v'(t) which brings the EOM into the form 


d (wa a 0 e ot cog Nt a 
wer ( a! ) — [x + W(t) a 90 ( —et2iwot cos Ot 0 )| ( o! ) : (8) 


So far the treatment of the parametric oscillator EOM has been exact. However, the present form of the equation 
allows to easily identify the terms responsible for parametric resonance. To this end, we apply the rotating wave 
approximation (RWA) (i) keeping in mind that the time-average of W(t) vanishes identically, and (ii) dropping any 
counter-rotating terms. Thus, we find that the dominant contribution to the dynamics appears for 2w9 = Q., which 
sets the critical driving frequency on resonance. The resulting effective RWA-EOM assumes the simple form: 


_d a’ _ Wo a a 
a (o) = (2 2) (). © 


Diagonalising the matrix on the right-hand side, we find the two Lyapunov exponents A1,2 = wo +iwoa/4. Hence, the 
maximum instability growth rate is set by awo/4. The stability criterion and the instability growth rate on resonance 
derived above by means of the RWA agree precisely with the standard results obtained using two-time perturbation 
theory or by other means [77]. 

In the following, we apply the same method and show explicitly the stability analysis for the driven BHM leading 
to the resonance condition Q. = 2Eaye(k) of the main text. Had it not been for the subindex aye, this result 
follows immediately from the above arguments for the parametric oscillator. We work in Bogoliubov theory which, 
as we show, captures the onset of instability. We begin by writing the rotating frame Hamiltonian of Eq. in 
momentum space and apply to it the Bogoliubov approximation. Parametrising the bosonic annihilation operator by 
by (t) = ux (t)b_(t = 0) — vt, (t)b"  (t = 0), with u(t = 0) = 1 and u(t = 0) = 0, the Heisenberg equations of motion 


read 
wa ue \) — [ e(k,t)+z2Jo+ Uno Uno Uk 
dt \ vu, J —Uno —le(—k, t) + zJp + Uno] Uk 


= Ga a +Uno esi ns Titittnes ) & ) 4 au s on ) (" ) - (10) 











where in the second line we separated the time average. Here no is the condensate fraction [of the time-averaged 
Hamiltonian], and the periodic function g,(t) = e(k,t) — €ave(k). It will prove convenient to first perform a static 
Bogoliubov transformation M;,(0), which diagonalises the time-averaged Hamiltonian: 


Uk \ u, \ _ { cosh(@,) sinh(@x) Uy 
( Uk ) Atel?) ( Vk ) 7 ( sinh(6x) cosh(Ox) ) ( Vk ) , ce 
The Bogoliubov angle is defined via cosh(20x) = (Eave(k) + n9U) /Eave(k) and sinh(20,) = n9U/Eave(k) with Eave(k) 
the corresponding Bogoliubov dispersion (see main text). The pseudounitary operator M,(@) has the property 
Mj (0)o7 My (0) = 07 {note that M-1(0) # MiL(8)]. 
The idea behind performing this Bogoliubov transformation is to bring the time-averaged Hamiltonian in diagonal 
form. At high-frequencies, the former represents the leading-order Floquet Hamiltonian [when expanded in powers 


of the inverse frequency], and thus this Bogoliubov transformation brings the state at time t = 0 to a basis which is 
close but not equal to the exact Floquet basis [finite Q~!-corrections are missing]: 


a fat \ z({ Uy, 
tae (ah) = Boer (1) 





: ( g(t) cosh? (Ox) + g(t) sinh? (A) 0 ) ( it, ) 
0 —g_(t) cosh? (Ox) — gx(t) sinh? (Ox) Uk 
+5 lout) +o-x(o]sinn(2m) (°, 4) (%). (12) 
Introducing the short-hand notation 
— ( gx(t) cosh? (x) + g_x(t) sinh? (Ox) 0 
W(t) ~ ( 0 —g_x(t) cosh? (Ax) = g«(t) sinh? (Ax) ) : 
h(t) = 5 lou(t) + 9-4 (0) (13) 
the EOM readily assumes the form: 
iS (UE) = enelloo® + Wiele] (Me) +sinnc2nyanto (9, 5) (GE). (14) 


Notice that the diagonal matrix W;,(t) has a zero time average, a property inherited from the function g(t). In the 
non-interacting limit, where one can integrate the EOM exactly, we have 0, — 0 and the time-dependent term W;(t) 
results in a trivial dynamical phase whose origin can be traced back to the kick operator of Floquet theory{5} [6]. 
Thus, we find that at small U, the kick operator attains additional contributions due to the Bogoliubov dressing in 
the diagonal matrix W(t). Although parametric instability is a phenomenon believed to originate in the micromotion 
itself [the stroboscopic Floquet Hamiltonian, if it exists as a local operator, is a hermitian operator and therefore has 
real eigenvalues], in both the weakly-interacting and the non-interacting case this W,(t)-term has no significance 
for the onset of the parametric instability, as we show below. We stress that W;,(t) is not the only contribution to 
the micromotion; part of the latter is due to the function h,(t). Thus, special attention must be paid to the term 
proportional to hy(t) sinh(20,) ~ noUh,(t). In the instantaneous diagonal basis, this results to the coupling of the 
form hy (ty cy, and is responsible for the drive-assisted scattering of bosons out of the condensate mentioned in the 
main text. 

In order to make the nature of parametric resonance visible, we perform yet another unitary time-dependent 
transformation %j,(t) = e7#e(tul (t), a, (t) = vi (t). This transformation is needed to bring the relative dynamical 
phases of uw, and v, with energy Eaye(k) at the same footing. In some physical sense, the onset of the parametric 
instability for a quantum harmonic oscillator is due to an enhanced mismatch of the relative dynamical phases 
accumulated by the operators b(t) and b!(t) in the presence of the periodic drive. The EOM now takes the form 


aft) — ; 0 hy(t)e~2#Eave(kt ii, 
“at ( ) - [Bat as W(t) 7 sinh (20x) ( —hy (t)e2#Bave(k)t 0 a, . (15) 


So far the above analysis was rather involved but exact. We can make further progress by Flourier-expanding 
the periodic function gx(t) = Jo Xizo ch. (Cje" [note that gx(t) has zero mean by definition] which leads to hy(t) = 


Jo/2 io [ch (C) +e, (¢)Je"*.. The Fourier coefficients cl, (¢) are closely related to Bessel functions. Now we apply the 
rotating wave approximation (RWA) to the above equation. This brings about a pronounced dominant contribution 
from the drive for Q, = 2Eayve(k) coming from the slowly rotating oscillatory off-diagonal terms. This behaviour is in 
precise agreement with the parametric resonance condition from the main text which was verified numerically using 
the WCCA. Recalling that the time-average of g(t) vanishes, the dynamics of the system on resonance is governed 
by the following set of effective equations 


Ls E,yo(k) Jo/2sinh(2.) [ee (¢) + =, (6)] 
= (at ) @retcs [eb (¢) + ch,(C)] Byye(k) ) (6) 


Diagonalising this matrix, we find the Lyapunov exponents on resonance [we avoid the terminology quasienergies since 
unlike quasienergies the Lyapunov exponents can be complex numbers], 








c1.2 = Fave(k) £4 sinh(26.)/ Ta + KO] EO + ARO], (a7) 


whose imaginary part is responsible for the exponential growth of the parametrically-unstable solution within Bogoli- 
ubov theory. Notice that sinh(20,) ~ noU defines the instability growth rate which is directly tied to the heating rate 
at short times in the parametrically unstable regime. At small driving amplitudes, the imaginary part of the Lyapunov 
exponent scales linearly with the drive strength ¢ = A/Q, as expected. Higher-order photon absorption resonances 
can be taken into account by using the higher-order Fourier coefficients e(¢ ). The behaviour of the system in the 
vicinity of the parametric resonance, on the other hand, can be analysed by introducing a small detuning 6 = Q—2F,. 





Finally, we remark that, the above analysis within Bogoliubov theory is not expected to produce the correct 
dynamics in the unstable regimes due to the lack of particle number conservation. Instead, one needs to further 
develop the perturbation theory by extending it to the WCCA. 
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FIG. 4: Comparison between the theoretical Bogoliubov stability boundary (black solid line) and the numerical solution to 
the WCCA equations discussed in the main text. Similarly to Fig. 1 of the main text, red circles mark the unstable while blue 
stars - the stable phase. The inset zooms in the small-¢ region to provide a better resolution. 


Figure |4|] shows a direct comparison between the stability criterion for the transition boundary 2 > 2Eaye(k) © 
“stable”, shown as a solid black line, and the numerical solution to the WCCA equations. Numerically, the condensate 
is defined to be stable whenever, after 800 periods of stroboscopic evolution, no decay is visible. In general, we find an 
excellent agreement between the two approaches; the only difference comes for the points lying right on the transition 
boundary which appear to be stable for small ¢ and unstable for large ¢. We attribute this to the uncertainty in 
determining the numerical phase boundary: indeed, coming from the stable phase the lifetime of the condensate can 
be very long, thus exceeding the 800 periods of evolution time. 
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DERIVATION OF THE EQUATIONS OF MOTION WITHIN THE WEAK-COUPLING CONSERVING 
APPROXIMATION 


In this appendix, we derive the equations of motion (EOM) using the weak-coupling conserving approximation 
(WCCA). We are interested in studying the periodically driven Bose-Hubbard model on a 2D lattice: 


H(t) =_ Dd Ji; (t)b'b; +h.c. + 5 Danis = 1), (18) 


In order to treat the spontaneous symmetry breaking of the condensate efficiently, we introduce the Bogoliubov 
spinor for the bosonic fields b > ba, with a = 1,2, where b; = b and by = b*. Adopting the notation (j,t) = x, the 
time-dependent action can be cast into the compact form 


Sb, b*] = So+ Sint 
Sob] = 5 f aari(a) (GLP) , vb) 
Sint [b, O°] = — yf aeedsc(e— yb" (28 (e)6(a)b(a) (19) 


where the integral over time is taken along the Keldysh roundtrip contour C and we introduced the delta 
function 6¢(# — x’) = d¢(t — t’)d;;". In Bogoliubov space, the noninteracting Green’s function has thus the form 


-1 10; + Jiz(t) 0 
(he), ( 0 18, + IE (t) ).: ey) 
We define the vacuum expectation value (VEV) v(x) and the quasiparticle (phonon) propagator G(a, y) as 
a(t) = (bCe))» 1Gas(esa) = (ba2)05(WN)e = ( OME. Cedelane (21) 


The microscopically occupied fields are denoted with a tilde b(x). Hence the Green’s function G defined above does 
not include the condensate fraction. The effective action is given by the double Legendre transform of the original 
action w.r.t. the VEV y() and the correlator Gap(x, y) [45) [79]: 


Tly,G] = Slp.g"] + 5TrllogG-] + 5TIGs"(y)G] - Tay, GI 
Slee] = f andy" @)G Lavo) -F fazio@l= fara @er(exvow), (22) 


where the sum over the Bogoliubov-Nambu index a is implicit. The Bogoliubov propagator Go ‘(x, Y;~) generates 
the motion of the Gross-Pitaevskii equation. Notice that it depends on the field y itself since the GPE is nonlinear. 
From that we obtain the inverse Bogoliubov propagator (Go")_, (#, yy) via: 


L pd fy Slee] 1 payed oO eee. 2p(x)|? p(x)? 
3 (Go = (x,y; Y) = dpe (x) dpy(y) ~ 5 (Gales (x,y) 2 da( y) ( (y(a)*)? Q(x)? }3 





(23) 


So far the calculation is exact, however, we have not specified the Luttinger-Ward functional ['2[y,G] yet which 
is the sum of all two-particle irreducible diagrams and thus has to be treated approximately. Here, we consider a 
weak-coupling expansion which amounts to consider diagrams to first order in U, see Fig. 

The EOM for the VEV and the propagator are obtained by making the effective action [ stationary with respect 


to the fields, oTly.C] = 0 and oT lei6 = 0, which lead to: 
7) Gab 





[ew (Gace) 11 (2, 9) 9(y) — Ue" (2) yp? (x) — U 29(2)Gu(z, 2) + y*(2)Gia(z,2)) = 0, 
iGy2 + gy? 


4 2(iGi1 + |y|?) Ny _+# 
» cor (t) — Udc(x — y) ( iGo + (y*)? 206Goo + lvl?) ) | Gocl(t,t) = dacdc(t—t), (24) 
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nes zc 
2x “15 Xv Gio(z, x) eo @ Goi (x, x) 
U 


=] 


Go(a, x) 2 


iF (0? x2 f dtoGrs(e,2)Goal0,2) -i5)? x | dteGre(e,2)Goi(e,2) 


FIG. 5: All two-particle irreducible diagrams which enter Iz to first order in U, with their proper combinatorial factors. The 
diagrams can be turned into equations using the following Feynman rules: (i) a factor of —iU/2 for each vertex, and (ii) a 
factor of i for each closed loop. By symmetry G11 = G22 and Giz = (Ga1)". 


and the Green’s function multiplication in the second equation above is understood in the matrix-multiplication sense: 
(AB)(a,z) = i. A(z, y)B(y, z). We remark that these EOM are equivalent to the Bogoliubov-Hartree-Fock EOM 


derived in Ref. [49] when starting from the lab frame Hamiltonian (see main text), making the ansatz b.=9 = y+ baa; 
bp40 = bezo, and then linearises any cubic terms in bx. 


Next, we open the closed time contour[79] by decomposing the Green’s function into a spectral part p(x, y) and a 
statistical part F(x, y) according to 


iG(a,2') = F(a,2') — 5P(@, 2" snc (t -?), 


rn _ 1 ny — 1 ( ({o(x),b'(a')})  ({b(a), &(a')}) 
Fase!) = (aCe) HG = 3 ( Ebay Blah) (BN Bee ), 
' ' b(x), bt (a! d(x), b(a’ 
pastor!) = iba) HheMDe= F(T Beryh (BHD Beh) = 
The following relations follow immediately from the above definitions: 
Fio(z,2') = Fio(z’, x), pie(x, 2’) = —pia(2’, 2), 
Foi (a,2') = Fai(2',x), poi(x, 2’) = —pai(2’, 2), 
Fyo(a, 2’) = Fey’), pi2(x, 2’) = pai (x, 2’), 
Fy (2, 2’) a2 Fi, (a, 2), pir(x, 2’) = —pj, (2’, 2), 
F22(x, 2’) = FY, (2, 2"), p22(x, 2’) = Pia (@,@ ). (26) 


We now assume that the system is translationally invariant, with periodic boundary conditions. We find the following 
system of coupled nonlinear EOM in momentum space for the condensate 


iOrp(t) = (zJo— pw) p(t) + ex=0(t) y(t) 


IZ lo)" [p@)l? + 2v(t) / Fu(t,t:@) + o@)" / Fralt,ta)| 


(27) 
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and the statistical correlator F’ 
i0;F 1 (t, t'; k) = oe — Lb) Fir( (t, t; ;k) +ex(t )Fii(t,t'; k) 
+7 [2(lewr + / Fultstea)) Falttsk) + ([oOF + f Falt,ta)) Fat esnyl], 
q 


10; Fyo(t, t'; k) = il (t, t; ;k) )+ez(t )\Fio(t,v’; k) 


nl? lo(t))? + [Pate ()) Faaltst k) + + (loo) P+ f Faltt a) [Fis (t, t’; k)] ik 


For completeness, we also give the equations of motion for the spectral correlators p which, on the other hand, obey 


(28) 


WOrpiu(t,tUsk) = (2Jo— pw)pult, tsk) + ex(t)pu(t, 5 k) 


+e [2 (Isto + / Fu(t,tsa)) ait Ds Co + Faa(tstsa)) laalttA)1 | 
Orpr2(t, Usk) = (zJo— p)pra(t, tsk) + ex(t)pra(t, 5 k) 


+o 2 Cai + [ Fulsta) pia(t, t's k) + Co + | Fiattta) lout 
(29) 


In the above equations, z is the coordination number, ¢,(t) is the time-periodic free dispersion in the rotating frame, 
and the integrals are all taken over the Brillouin zone. 

Furthermore, if one is interested in the equal-time correlation of the statistical correlator F’, using the symmetry 
relations in Eq. one arrives at the somewhat simplified equations 


Fiu(t,tk) = zim { 5 (iv + [ Palit a) [Fio(t, t; k)]* 1, 


i0:Fya(t,t;k) = 24 (24 — p)Fio(t, t;k) + en (t) Fie(t,t; k) 
+ [2(loor+ / Fu(tstsa)) Fralttsk) + ((eOP + / Faltytea)) Fatty] }. (30) 


WCCA Equations of Motion for the BHM on a General Bipartite Lattice 


We now generalise the WCCA EOM to any bipartite lattice. Consider a bipartite lattice with the two sublattices 
labelled by A and B and periodic boundary conditions. Each sublattice contains N4 = Ng = N,/2 number of sites. 
With this definition, additionally to the Bogoliubov index a = 1, 2, all correlators carry an additional index a = A, B, 
and so does the condensate fraction. 

The extended system of equal-time equations of motion for the WCCA of the BHM reads 


ip (t) = (zJo— nw) pA (t) + ex=o(t)e? (t) 
oe [aah (oP +264) [ rAAeng + loroy" f FiAtto) | 


Na 
iOp? (t) = (zJo— we? (t) + Exao(t)y4 (t) 
U By4)]* [, ,Brp\]2 B BB By4y)* BB 
+e [Pol lO) +2080 fre wo+ oO)" [FBP eO), (31) 


AFA (tthk) = atm {ext FAP em) +e (eco) + / Fito) nanny}, 


A 


FRM tk) = Bn eeORAM EAR) + ae (LOMO) + f PEM) [REMC 
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U 
iOFiV? (t,t)k) = ex(t) (FAP (66%) — PAA G&)) + ml 


+2 (Ie? — |p? (t)P? + / Fi (t, tq) - FA (t,150)) Fe GER) 
qd 


_ ([e*o]? +/ FiAt.ta)) [FAP (t,t:k))" - ([67(0)]” i: / FBP (ta) FAP (tt i). 


q 





iF AA (t,t34) = 24 (24 — p)FAYA (t, t3k) + ee FIBA (tk) 

ty 2(leMOr+ [ rieeo) Avene + ([otol+ f ratGen) Feeney |}, 
iF BP (t,t:k) = 24 ee(OFAP Ct k) + (2Jo — WFR" (tts k) 

ty, P(e + [Pa eeo) HERG + (fer + [ ABPGe a) [Ren J, 
iF in? (t,t3k) = 2(2Jo — w)Fi2" (t,t; k) + (en (t)" Fin" (tts k) + ex (F127 (tt5h)) + ma | 


+2 (1s? + |p? (t)/? + ‘ Fah (t, tig) + FEF(.t4) ie Gk) 
qd 


r ([e*o]? +/ FiAt.ta)) [FAP (tk) - ([670)]” a i FBP (ta) FAP (tt i). (32) 


qd 





All integrals in nape and are taken over the reduced Brillouin zone (w.r.t. the AB-sublattice symmetry). 
Equations and (32) constitute a coupled set of non-linear equations, the solution of which produces the dynamics 
discussed in the main text. Note that these EOM can be applied to systems with arbitrary time-dependence (not 
necessarily a periodic one) and on an arbitrary bipartite lattice, such as the honeycomb lattice. 

For the analysis in the main text, the initial condition for the condensate fractions is chosen to be |y“(0)|?/Ns = no/2 
for ¢ < 2.405 and |y?(0)|?/N, = no/2 for ¢ > 2.405, where no is the total condensate fraction for the non-driven 
model in Bogoliubov theory. 


VALIDITY OF THE WCCA AND THERMALISATION TIMESCALES 


In this Appendix we estimate the timescales on which the WCCA gives a reliable description of the physics, and 
discuss the dominant processes that (in the weak coupling regime) destabilize the various pre-thermal steady states 
discussed in the main text. We discuss each of the three regimes separately. 

Parametrically unstable region. In this regime, the prethermalized phase is the one in which the momentum 
distribution is sharply peaked along momentum-space arcs as in Fig. 3 (a) (see main text). As in the main text, 
we treat the time-averaged dispersion as the unperturbed Hamiltonian and look at the parametric instability growth 
rate growth rate due to a perturbation of the form JoTi(C) bi bee. The matrix element for pair creation is then 
~ UnoJoJi(C)/Wave (the Fourier coefficient c!(¢) from the parametric instability analysis above is essentially given 
by the Bessel function), and for reasonably large drives this is linear in U. Parametric instability predicts that these 
features will grow at the rate [ ~ Ung, where no is the condensate amplitude. The decay rate (i.e., inverse lifetime) 
of the quasiparticles along these arcs, once they are formed, is limited by collisions, and Fermi’s Golden Rule implies 
that this decay rate is of order U?; this is the rate at which these features spread out in momentum space. Thus there 
is a parametric separation in U between the formation and decay rate of these peaks. The leading collisional process 
comes from cubic terms of the form Up bh. bis by, -k (plus appropriate conjugates) in the Hamiltonian. The Golden 
Rule rate for this particular process is 


T.(k) ~ U? non N 2p (Eave(k)) (33) 


where N2p(Eave(k)) ~ f dE’d?q5(E’ — Eave(q))5(Eave(k) — E’ — Eave(q —k)) is the accessible two-particle density 
of states. Here, Eaye(q) is the energy of an excitation with quasimomentum q. On dimensional grounds this two- 
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particle density of states must be inversely proportional to Waye; thus, the overall Golden Rule lifetime of a particular 
quasiparticle state will go as 


I. (k) rm U? nonk/Wave (34) 


up to a multiplicative constant. The ratio between the decay rate and the creation rate scales as Un, / (Joi (C)). Thus 
the decay rate of a mode is slower than the creation rate whenever the condensate amplitude is large compared with the 
population of the mode (essentially because the matrix element is not Bose-enhanced to the same degree). However, 
the decay rate is also suppressed with decreasing the interactions (expected) or increasing the drive amplitude. At 
short times, when the condensate is not appreciably depleted, the WCCA is therefore reliable; however, when the 
depletion becomes large the WCCA also fails. Thus the regimes of validity of the WCCA and Bogoliubov theory in 
the parametrically unstable regime are essentially the same, although the WCCA has the advantage of respecting 
particle number conservation exactly at all times. 


Stable region. In the stable region there are two types of physical processes beyond the WCCA. (1) The excitations 
created by the original quench into the phase have finite collisional lifetimes, as discussed above. The momentum- 
space patterns in the stable region will dephase on this Golden-Rule timescale [.; however, the condensate fraction 
will remain large and stable even after dephasing. (2) Eventually, the system will absorb energy from the drive. If the 
drive frequency is 2“ and the bandwidth of single-particle excitations is Waye, then resonant absorption must involve 
at least m = Q/Waye quasiparticles. It is straightforward to check that the associated Golden Rule rate, at weak 
coupling, is of the form U™/W7>'. When U is sufficiently small, this heating timescale is much longer than the 
timescale on which the momentum-space patterns dephase; thus the system should remain stable for extremely long 
times at high frequencies. 


Dynamically unstable region. In this regime, the growth rates of unstable modes are of order \/ Jaye(¢)U, whereas 
the collision rates are of order U? /Wave at best, so at weak coupling we have a parametric window in U where the 
WCCA remains valid. 


PHASE TRANSITION REGION AROUND ¢ = 2.405 


In this appendix, we discuss the dynamics governed by the WCCA close to the first zero of the Bessel function, 
¢ = 2.405, where the dispersion of the 2: > co Hamiltonian becomes flat (central grey region in Fig. 1, main text). 
For ¢ < 2.405 the dispersion of the free theory U = 0 supports a stable minimum for k = 0, while for ¢ > 2.405 the 
stable minimum appears at k = 7. Since the the two stable regions support different momentum modes, a phase 
transition occurs in between them. Therefore, it is required that one allows for a macroscopic population of both the 
modes in the immediate vicinity of ¢ = 2.405. 


This can be achieved by reducing the translational symmetry of the problem. Intuitively, a condensate at k = 7 
with amplitude y,—, flips a sign on every other site. Hence, one can choose to work in the original (momentum- 
resolved) basis (~p=0,r—x), Or in the site-resolved basis (y4,y?). The two are related by a rotation. In order 
allow for a dynamical population of the yz=, condensate, Eqs. (31) and require that the initial condition for 
yr (0) = 1/V2 (~4(0) — y?(0)) be nonzero. In the AB-basis, this is equivalent to saying that there is a slight difference 
in the condensate occupation on the two sublattices. Physically, this imbalance is caused by spontaneous symmetry 
breaking. However, in the WCCA one has to put in this imbalance by hand. In the following we refer to the small 
value s = |y7(0)|? as seed. 


When the effective dispersion becomes flat ¢ ~ 2.405 [Fig [6], the condensate undergoes oscillations between the 0 
and w modes, with a period ~ 1/U for small U. This behaviour is reminiscent of the collapse-and-revival effects seen 
for a BEC that is suddenly quenched into the Mott insulating phase[80}, although the dynamics governed by WCCA 
is classical. The period of the transfer oscillations is also seed-dependent and increases with s > 0. Even though our 
approximation does not capture a true Mott insulating phase, the nonlinearities included in the WCCA are sufficient 
to give rise to these oscillations. Physically speaking, a quasiparticle-mediated channel is opened, through which 
particles flow from the condensate at k = 0 tok = z. Although it is present at any ¢, this channel is only effective 
when the dispersion is sufficiently flat since the amplitude for the phonon- mediated transition yroo > bf — YPkren 
scales as (U/Jo)°. 
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FIG. 6: Time evolution of the condensate fractions for ¢ = 2.405 and Q/Jo = 20 starting from a Bogoliubov initial state 
localised at k = 0. The seed size is s = 1% and U/Jo = 0.2. 
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FIG. 7: Total energy density of condensate and quasiparticles as a function of time for U/Jo = 0.2 following the quench with 
frequency and amplitude as stated in the legend. 


TIME-DEPENDENCE OF THE ENERGY 


Last, we briefly address the issue of heating. Fig. [7|shows the excess total (i.e., condensate plus quasiparticle) energy 
density in the system, relative to the non-driven state. Due to the abrupt turn-on of the periodic circularly polarised 
modulation, the energy changes discontinuously at t = 0. As expected, the energy density increases due to heating 
in the parametrically unstable region, saturates in the stable region, and exhibits a small growth for ¢ = 3.8. Notice 
the different behaviour in the parametrically unstable region compared to the dynamically unstable one: while in the 
former the energy grows due to the population of modes lying on the high-energy surface, in the latter the dynamically 
unstable modes appear close enough to the origin [cf. Fig. 3, panel (c) in main text] so that the growth in energy density 
past the quench value is not substantial. Note that the system does not heat up even at fairly long times whenever 
the parameters are chosen to be in the stable region of the stability diagram. Although ergodic periodically-driven 
systems are expected to eventually heat up to infinite temperature|33} [38], in the weak-coupling limit this 
heating timescale (which is due to collisions between quasiparticles) is parametrically slower in the “stable” regimes 
of our phase diagram than in the “unstable” regimes. Thus, for a range of present-day experiments, we expect that 
in the stable high-frequency regime there is no significant heating on experimentally relevant timescales. 
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TIME-DEPENDENCE OF THE MOMENTUM DISTRIBUTION FUNCTION 


For the time-evolution of the momentum distribution function, we refer to the three videos in the supplementary 
material. The lower left panel shows the quasiparticle momentum distribution over the first Brillouin zone, while 
the upper left panel is a top view of the same. The upper right panel displays the time-evolution of the condensate 
fraction, while the lower right panel shows the energy density. The parameters for each simulation can be found in 
the title. The three movies correspond to the points in the stability diagram marked by (a), (b), (c) in Fig. 1 (main 
text). 





